Classify LAS.sml

##Thin LAS File.sml
#  This script is designed to reclassify LiDAR point clouds to Ground and non-ground based on distance from a DEM surface.
#
#  Script reads LAS records in and determines location of point on DEM raster. If z value of point is different by more than 
#  threshold the point is classified as "Unclassified" if the values is close to the DEM value then the point is classified as "Ground".
numeric threshold = 0.1;
# 
# Terranean Mapping Technologies, Brisbane, Australia: 
# 21/5/2009
#####################################################################
clear();
class RVC_OBJITEM InputLASList[];
class STRINGLIST listkeys;
class RVC_RASTER DEM;
class RECT3D extents;
class RVC_GEOREFERENCE georef;
class SR_COORDREFSYS refsys;
class TRANS2D_MAPGEN demMapToObject;
class RVC_SHAPE lasshape;
class RVC_DESCRIPTOR descriptor;
class RVC_DBASE_SHAPE lasDB;
class RVC_DBTABLE lasTable;
class RVC_DBTABLE_RECORD lasRecord;
class POINT2D lasMapPoint, demObjPoint;
string str$;
numeric l, nrecords, i, c, z;
print("Getting Input/Output datasets");
#Select Las files to process
DlgGetObjects("Select LAS File(s)", "ALL", InputLASList, "ExistingOnly");
listkeys = InputLASList.GetKeys();
#Select input filtered DEM
GetInputRaster(DEM);
DEM.GetExtents(extents);
georef.OpenLastUsed(DEM);
refsys = georef.GetCoordRefSys();
#Set up MaptoObject Transform
georef.GetTransparm(demMapToObject, 1, georef.GetCalibModel());
demMapToObject.InputCoordRefSys = refsys;
#Loop through all the selected LAS files and reclassify ground points based threshold distance from DEM surface
print("Classifying LAS files");
for l = 1 to listkeys.GetNumItems() {
	#Open up the current LAS file and its main table.
	lasshape.Open(InputLASList[l], "Read");
	descriptor = lasshape.GetDescriptor();
	printf("\t%d - %s\n",l, FileNameGetName(descriptor.GetFullName()));
	lasDB.OpenAsSubobject(lasshape, "Write");
	lasTable.Open(lasDB, descriptor, "Write");
	#Loop through records and change classification to unclassified(1) for non ground points and ground(2) for ground points within a threshold of the DEM surface
	nrecords = lasTable.GetNumRecords() - 1;
	for i = 0 to nrecords{
		if(i % 10000 ==0){
			str$ = sprintf("Point %d of %d", i, nrecords);
			SetStatusMessage(str$);
			SetStatusBar(i, nrecords);
		}
		lasTable.Read(i, lasRecord);
		#read. x, y, z, and classification value from the Las database
		c = lasRecord.GetValue("Classification");
		lasMapPoint.x = lasRecord.GetValue("X");
		lasMapPoint.y = lasRecord.GetValue("Y");
		z = lasRecord.GetValue("Z");
		demObjPoint = TransPoint2D(lasMapPoint, demMapToObject);
		if(abs(DEM[demObjPoint.y+1, demObjPoint.x+1] - z) > threshold) c = 1;
		else c = 2;
		lasRecord.SetValue("Classification", c);
		lasTable.Write(i, lasRecord);
		nextrecord:
	}
	}
theend:
print("Done!");
#eof